This page contains the code to make antigenic maps and antibody landscapes. You can open the R Markdown file in RStudio to have an interactive step-by-step document and chose the “Visual” display instead of the “Source” display if you prefer it.
The data used here was published in Roessler et al., Nat Comm (2022), publicly available in the paper’s GitHub repository and in this GitHub repo. If you use the data, please cite the original publication. This data contains antibody neutralization titers. However, maps can also be made from binding data or hemagglutinin inhibition data.
If you use materials from this page, please use this repository’s DOI as reference and cite the Racmacs and/or ablandscapes package.
Scientific background on antigenic cartography can be found in Smith et al., Science (2004) and on antibody landscapes in Fonville et al., Science (2014).
A detailed introduction to the R package for making antigenic maps is given on the Racmacs reference page.
To run the Racmacs package and construct antigenic maps and antibody landscapes, we use R: R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
For examples of map diagnostic steps you can look at the Supplementary Material in Roessler et al., Nat Comm (2022) and Roessler et al., Nat Comm (2023).
Antigenic maps were constructed using the Racmacs package, Version 1.1.35: Wilks S (2022). Racmacs: R Antigenic Cartography Macros. https://acorg.github.io/Racmacs, https://github.com/acorg/Racmacs.
Antibody landscapes were constructed using the ablandscapes package, Version 1.1.0: Wilks S (2021). ablandscapes: Making Antibody landscapes Using R. R package version 1.1.0, https://github.com/acorg/ablandscapes.
We recommend calculating geometric mean titers and fold changes using the titertools package to deal with censored data: Wilks SH (2022). titertools: A statistical toolkit for the annalysis of censored titration data. R package version 0.0.0.9001.
If the slideshow does not render for you, try opening it in another
browser. It should open in Chrome. Alternatively, you can open the pdf
./Introduction_to_antigenic_cartography.pdf. This slideshow
gives a brief overview of what is demonstrated with real world data in
the remaining document.
# these are the packages used in the notebook and should be installed for it to run smoothly
library(tidyverse)
library(Racmacs)
library(ablandscapes)
library(Racmacs)
library(dplyr)
library(r3js)
library(knitr)
library(kableExtra)
library(patchwork)
#install titertools from GitHub: https://github.com/shwilks/titertools.git if you want to estimate <LOD values for GMT calculation
# library(titertools)
Antigenic maps should be created from single exposure sera or baseline vaccination sera. With just this type of sera (e.g. only 2 dose vaccinated), the map will be poorly triangulated, meaning the positional resolution of variants will be very low. Sera will not be positioned accurately, or at all if titers against the majority of variants are below the limit of detection (<LOD). Homologous sera from the different variants are needed to add resolution. The map by Roessler et al. (Nat Comm 2022) contains many different human first infection sera.
Importantly, antigenic maps always reflect the data that was used to construct them. So two maps with different types of input sera and variants may look different. However, if enough sera and variants are present in the input data to triangulate the items in euclidean space, maps from different source data should be consistent. Hence maps that look different might be an accurate representation of the input data, but not an accurate representation of antigenic space due to lack of diversity in the input (demonstrated below). To assess the suitability of input data for antigenic maps, map diagnostic steps are important.
Set the working directory to the root directory of your project for convenience.
Antigenic maps deal with titers below the limit of detection (<LOD) by allowing the euclidean map distance to be larger than or equal to the titer target distance to the LOD. This utilizes the information that is contained in this measurement: We know the distance is at least as much as the maximum titer to LOD distance. Here, the LOD was 16.
The input for making a map is a table of measured titers. The rows correspond to the titrated antigens, the columns to the measured sera. Here, I will extract a table from an existing map and demonstrate how to create a map from this table.
working_dir <- "~/Documents/smith/labbook/antigenic_cartography_introduction/"
# set working directory
setwd(working_dir)
# read in data for multiple exposure landscapes
multi_exposure_data_b <- read.csv("data/titer_data/titer_data_long.csv") %>%
mutate(sr_group = gsub("\\/alpha\\+E484K", "", sr_group))
# go to Additional section to find out more
multi_exposure_data <- multi_exposure_data_b %>%
filter(ag_name != "P.1.1")
# Racmacs function to read in antigenic map of .ace format
og_map <- read.acmap("data/maps/map-OmicronI+II+III-thresholded-single_exposure-P1m1.ace")
# extract Titer Table from map
titer_data_og <- titerTable(removeAntigens(og_map, "P.1.1"))
titer_data <- apply(titer_data_og, 2, function(x){
temp <- as.character(round(as.numeric(x)))
temp[is.na(temp)] <- "<16"
return(temp)
})
rownames(titer_data) <- rownames(titer_data_og)
# shorten colnames of samples
colnames(titer_data) <- sapply(colnames(titer_data), function(x){
paste(str_split(x, "_")[[1]][1:2], collapse = "_")
})
colnames(titer_data) <- gsub("\\/alpha\\+E484K", "", colnames(titer_data))
# read in colors for map
sr_group_colors <- read.csv("data/metadata/sr_group_colors.csv", sep = ";")
mapColors <- read.csv(file = './data/metadata/map-colors.csv', row.names = 'Antigen', header = TRUE)
mapColors[sr_group_colors$SerumGroup, "Color"] <- sr_group_colors$Color
# show titer table
kable(titer_data, "html", align = "r") %>%
kable_styling("striped", full_width = F)
| G21_delta conv. | G22_delta conv. | G23_delta conv. | G24_delta conv. | G25_delta conv. | G26_delta conv. | G27_delta conv. | F620_alpha conv. | F628_alpha conv. | F633_alpha conv. | F635_alpha conv. | F647_alpha conv. | F650_alpha conv. | F658_alpha conv. | F661_alpha conv. | F662_alpha conv. | F667_alpha conv. | C701_beta conv. | C709_beta conv. | C711_beta conv. | C715_beta conv. | C770_beta conv. | C850_beta conv. | C859_beta conv. | C860_beta conv. | G30_mRNA1273/mRNA1273 | G33_mRNA1273/mRNA1273 | G34_mRNA1273/mRNA1273 | G36_mRNA1273/mRNA1273 | G37_mRNA1273/mRNA1273 | G38_mRNA1273/mRNA1273 | G39_mRNA1273/mRNA1273 | G40_mRNA1273/mRNA1273 | G41_mRNA1273/mRNA1273 | G42_mRNA1273/mRNA1273 | E916_AZ/AZ | E918_AZ/AZ | E993_AZ/AZ | E995_AZ/AZ | E998_AZ/AZ | E999_AZ/AZ | E1000_AZ/AZ | F5_AZ/AZ | F6_AZ/AZ | F9_AZ/AZ | E919_AZ/BNT | E920_AZ/BNT | E921_AZ/BNT | F41_AZ/BNT | E994_AZ/BNT | E996_AZ/BNT | E997_AZ/BNT | F1_AZ/BNT | F2_AZ/BNT | F3_AZ/BNT | F110_BNT/BNT | F120_BNT/BNT | F121_BNT/BNT | F122_BNT/BNT | F124_BNT/BNT | F262_BNT/BNT | F263_BNT/BNT | F269_BNT/BNT | F271_BNT/BNT | F273_BNT/BNT | F289_BNT/BNT | G291_BA.1 conv. | G347_BA.1 conv. | G348_BA.1 conv. | G353_BA.1 conv. | G354_BA.1 conv. | G355_BA.1 conv. | G356_BA.1 conv. | G359_BA.1 conv. | G360_BA.1 conv. | G362_BA.1 conv. | G363_BA.1 conv. | G366_BA.1 conv. | G369_BA.1 conv. | G371_BA.1 conv. | G381_BA.1 conv. | G393_BA.1 conv. | G647_BA.1 conv. | G650_BA.1 conv. | G665_BA.2 conv. | G668_BA.2 conv. | G669_BA.2 conv. | G670_BA.2 conv. | G671_BA.2 conv. | G672_BA.2 conv. | G673_BA.2 conv. | G674_BA.2 conv. | G751_BA.2 conv. | G776_BA.2 conv. | G780_BA.2 conv. | G782_BA.2 conv. | 218 (Ischgl 1)_WT conv. | 220 (Ischgl 1)_WT conv. | 224 (Ischgl 1)_WT conv. | 260 (Ischgl 1)_WT conv. | 262 (Ischgl 1)_WT conv. | 278 (Ischgl 1)_WT conv. | 279 (Ischgl 1)_WT conv. | 280 (ischgl 1)_WT conv. | 298 (ischgl 1)_WT conv. | 299 (Ischgl 1)_WT conv. | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| D614G | 559 | 901 | 85 | 114 | <16 | <16 | 43 | 255 | 27 | 56 | 98 | 21 | 23 | 19 | 148 | 171 | 120 | 307 | <16 | <16 | 22 | 139 | 150 | 529 | <16 | 926 | 442 | 1631 | 1518 | 589 | 180 | 923 | 1185 | 97 | 348 | 198 | 36 | 325 | 649 | 52 | 43 | 140 | 120 | 839 | 148 | 1186 | 859 | 2106 | 915 | 2403 | 1003 | 512 | 1244 | 3582 | 1594 | 1159 | 1204 | 1495 | 654 | 1988 | 357 | 1592 | 560 | 2567 | 981 | 689 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | 17 | <16 | 20 | <16 | <16 | 22 | <16 | <16 | <16 | <16 | <16 | <16 | 26 | 788 | 204 | 17 | 165 | 200 | 463 | 123 | 132 | 75 | 63 | 200 | 132 | 124 |
| B.1.1.7 | 415 | 665 | 66 | 153 | <16 | <16 | 41 | 353 | 167 | 100 | 479 | 38 | 158 | 184 | 533 | 551 | 432 | 341 | <16 | <16 | 48 | 642 | 37 | 173 | <16 | 867 | 136 | 398 | 470 | 363 | 41 | 389 | 856 | 80 | 425 | 149 | 36 | 197 | 218 | 58 | 63 | 97 | 79 | 948 | 137 | 1864 | 868 | 1823 | 404 | 1642 | 1446 | 593 | 514 | 3849 | 843 | 923 | 531 | 1723 | 369 | 1020 | 315 | 1738 | 793 | 601 | 361 | 246 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | 36 | <16 | <16 | <16 | <16 | 54 | <16 | <16 | <16 | <16 | 49 | <16 | <16 | 448 | 247 | <16 | 155 | 60 | 357 | 228 | 32 | 88 | 47 | 200 | 58 | 66 |
| B.1.1.7+E484K | 16 | 450 | 96 | 117 | <16 | <16 | 70 | 345 | 142 | 56 | 361 | <16 | 130 | 75 | 265 | 266 | 401 | 789 | 137 | <16 | 72 | 598 | 139 | 723 | <16 | 441 | 178 | 556 | 248 | 192 | 96 | 368 | 397 | 41 | 236 | 74 | 17 | 281 | 124 | 26 | 35 | 63 | 76 | 79 | 130 | 537 | 347 | 424 | 171 | 580 | 402 | 302 | 346 | 4433 | 719 | 1269 | 625 | 1272 | 213 | 796 | 353 | 754 | 81 | 473 | 279 | 240 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | 128 | <16 | 20 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | 40 | <16 | <16 | 280 | 210 | <16 | 283 | 38 | 320 | 392 | 23 | <16 | 36 | 26 | 54 | <16 |
| B.1.351 | 17 | 1396 | 164 | 33 | <16 | <16 | 40 | 156 | 68 | 54 | 429 | <16 | 219 | 137 | 114 | 205 | 305 | 909 | 102 | <16 | 533 | 2014 | 634 | 2025 | <16 | 130 | 104 | 614 | 227 | 75 | 24 | 225 | 330 | 45 | 279 | 147 | <16 | 251 | 300 | 71 | 30 | 47 | 87 | 67 | 91 | 512 | 457 | 466 | 199 | 709 | 387 | 477 | 226 | 796 | 397 | 235 | 108 | 643 | 778 | 325 | 133 | 615 | 68 | 185 | 176 | 116 | 27 | <16 | <16 | <16 | <16 | <16 | <16 | 24 | <16 | 25 | <16 | <16 | <16 | <16 | 25 | <16 | <16 | <16 | <16 | 39 | <16 | <16 | <16 | <16 | 47 | <16 | <16 | 178 | 97 | <16 | 228 | 108 | 178 | 87 | 22 | <16 | 22 | 25 | 45 | <16 |
| B.1.617.2 | 448 | 632 | 126 | 2726 | 227 | 91 | 146 | 78 | 50 | 77 | 89 | <16 | 18 | 17 | 43 | 116 | 133 | 39 | <16 | <16 | <16 | 111 | <16 | 254 | <16 | 713 | 67 | 484 | 414 | 243 | 31 | 400 | 314 | 24 | 107 | 181 | 25 | 122 | 274 | 191 | 24 | 44 | 101 | 182 | 129 | 716 | 639 | 873 | 291 | 895 | 536 | 294 | 185 | 864 | 458 | 237 | 178 | 349 | 182 | 344 | 155 | 293 | 191 | 384 | 195 | 126 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | 47 | 30 | <16 | <16 | 53 | <16 | <16 | 20 | <16 | 178 | <16 | <16 | 21 | <16 | 20 | 36 | <16 | 44 | <16 | 47 | 231 | 275 | <16 | 257 | 38 | 569 | 259 | 94 | 93 | 36 | 223 | 48 | 51 |
| BA.1 | <16 | <16 | <16 | 17 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | 26 | <16 | <16 | <16 | <16 | <16 | <16 | 16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | 62 | <16 | 32 | <16 | 42 | 25 | <16 | <16 | 69 | 23 | 28 | <16 | 19 | <16 | 50 | <16 | 65 | <16 | 18 | <16 | <16 | 106 | 31 | 68 | <16 | 44 | 21 | 34 | 122 | 76 | 54 | 144 | 42 | 51 | 80 | 140 | <16 | 365 | 52 | 30 | <16 | <16 | 21 | 28 | <16 | <16 | <16 | 45 | 111 | <16 | <16 | <16 | <16 | 109 | <16 | <16 | <16 | <16 | <16 | <16 | 22 |
| BA.2 | 30 | 113 | 58 | 34 | <16 | <16 | 83 | 24 | 17 | <16 | 33 | <16 | <16 | <16 | 18 | 56 | 98 | 24 | <16 | <16 | 30 | 100 | 20 | <16 | <16 | <16 | 85 | 273 | 112 | 184 | <16 | 115 | 137 | <16 | 120 | 46 | <16 | 34 | <16 | <16 | <16 | 147 | 48 | 30 | 38 | 257 | 146 | 351 | 162 | 440 | 125 | 217 | 1024 | 632 | 160 | 254 | 151 | 501 | 82 | 273 | <16 | 186 | <16 | 309 | 95 | 109 | 17 | <16 | <16 | <16 | <16 | <16 | 46 | 21 | 23 | 37 | 51 | 113 | 27 | 37 | 204 | <16 | 693 | 44 | 391 | 88 | 1101 | 553 | 335 | 96 | 157 | 120 | 525 | 548 | 69 | 1055 | 40 | 41 | 170 | 149 | <16 | <16 | <16 | <16 | 25 | <16 |
| BA.5 | <16 | 447 | 68 | 67 | <16 | <16 | 21 | 21 | <16 | <16 | 52 | <16 | <16 | <16 | 17 | 141 | 233 | <16 | <16 | <16 | 26 | 177 | <16 | <16 | <16 | 18 | 23 | 91 | 84 | 25 | <16 | 29 | 53 | <16 | 55 | <16 | <16 | <16 | <16 | <16 | <16 | 21 | <16 | <16 | 26 | 176 | 56 | 202 | 38 | 67 | 60 | 57 | 66 | 343 | 50 | 47 | 28 | 36 | 38 | 83 | 28 | 68 | <16 | 91 | 22 | 47 | <16 | <16 | <16 | <16 | <16 | <16 | <16 | 29 | <16 | <16 | <16 | <16 | <16 | 26 | 28 | <16 | 33 | <16 | 56 | 50 | 17 | 52 | 41 | 37 | 66 | 24 | 69 | 447 | 47 | 27 | 74 | <16 | 72 | 62 | <16 | <16 | <16 | <16 | 26 | <16 |
Let’s have a look at the titer data:
multi_exposure_data %>%
# log2 transform the data and set <LOD values to LOD/2
mutate(logtiter = ifelse(titer == "<16", log2(16/20), log2(as.numeric(titer)/10)),
sr_group = factor(sr_group, levels = sr_group_colors$SerumGroup)) -> data_long
# calculate
data_long %>%
group_by(
sr_group,
ag_name
) %>%
# if you want to use quick LOD/2 method for GMT calculation
summarize(logtiter = mean(logtiter, na.rm = TRUE)) %>%
# We recommend using the titertools package for GMT/fold change calculation where <LOD values are estimated in
# a bayesian approach. As it runs a bayesian model in the background, it takes a bit longer to run. The code to run
# it is commented out below. It returns the mean, lower and upper estimates of the mean and of the standard deviation
# summarize(gmt = titertools::gmt(titer, dilution_stepsize = 0)["mean", "estimate"]) %>%
# manually set GMT's that are lower than that to LOD2
mutate(logtiter = ifelse(logtiter < log2(0.8), log2(0.8), logtiter),
sr_name = "GMT",
gmt = logtiter)-> gmt_data
do_titerplot <- function(data_long, target_sr_groups, gmt_data, sr_group_colors, ag_order){
data_long %>%
filter(sr_group %in% target_sr_groups) %>%
mutate(logtiter = ifelse(titer == "<16", log2(16/20), log2(as.numeric(titer)/10))) %>%
ggplot(aes(x = ag_name, y = logtiter, color = sr_group, group = sr_name)) +
geom_line(aes(group = sr_name), alpha = 0.3) +
geom_point(alpha = 0.3) +
geom_line(data = gmt_data %>%
filter(sr_group %in% target_sr_groups), size = 1.5) +
geom_point(data = gmt_data %>%
filter(sr_group %in% target_sr_groups), fill = "white", shape = 21, size = 3) +
scale_color_manual(values = setNames(sr_group_colors$Color, sr_group_colors$SerumGroup),
name = "Serum group") +
scale_x_discrete(limits = ag_order,
name = "Variant") +
scale_y_continuous(name = "Titer",
labels = function(x) round(2^x*10),
breaks = seq(log2(16/10), 12, 1)) +
annotate("rect",
xmin = -Inf,
xmax = Inf,
ymin = -Inf,
ymax = log2(16/10),
fill = "grey",
alpha = 0.6,
color = NA) +
facet_wrap(~sr_group) +
theme_bw() +
theme(strip.background = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)) -> p
return(p)
}
single_exposure_groups <- c('delta conv.', 'alpha conv.', 'beta conv.', 'mRNA1273/mRNA1273', 'AZ/AZ', 'AZ/BNT', 'BNT/BNT', 'BA.1 conv.', 'BA.2 conv.', 'WT conv.')
# plot the data
do_titerplot(data_long, target_sr_groups = single_exposure_groups,
gmt_data = gmt_data,
sr_group_colors = sr_group_colors,
ag_order = rownames(titer_data))
Individual titers are shown as faint lines, geometric mean titers (GMT) as thick lines. We can see some trends, for example BA.1’s escape in most serum groups.
# set seed for random optimization start
set.seed(100)
# create the acmap object
map <- acmap(
ag_names = rownames(titer_data),
sr_names = colnames(titer_data),
titer_table = titer_data)
# The dilution stepsize gives the dilution steps of the neutralization assay on the log2 scale, e.g.: 0 = continuous read out, 1 = two-fold dilutions
dilutionStepsize(map) <- 0
# optimize the map in 2 dimensions
map <- optimizeMap(
map,
number_of_dimensions = 2,
number_of_optimizations = 1000,
minimum_column_basis = "none",
options = list(ignore_disconnected = TRUE) # Map contains disconnected points (points that are not connected through any path of detectable titers so cannot be coordinated relative to each other)
)
plot(map, plot_labels = "antigens")
Above is a basic map. We will add meta-information in the next step to make it interpretable:
# Add fill color for map
agFill(map) <- mapColors[agNames(map),]
agGroups(map) <- agNames(map)
# Add serum group = serum cohort to map. It is stored here in the sample name, the second field after "_"
sr_groups <- unlist(lapply(srNames(map), function(x) str_split(x, "_")[[1]][2]))
srGroups(map) <- factor(sr_groups, levels = sr_group_colors$SerumGroup)
srOutline(map) <- mapColors[as.character(srGroups(map)),]
# Add general style guidelines to map
srOutlineWidth(map) <- 1
srSize(map) <- 9
agSize(map) <- 18
srOutline(map) <- adjustcolor(srOutline(map), alpha.f = 0.7)
# make some antigens smaller than others
sublineages <- c("B.1.1.7+E484K")
agSize(map)[agNames(map) %in% sublineages] <- 15
# align map to example map. Important: antigen names have to match
map <- realignMap(map, og_map)
# get map limits for plotting
lims_no_zoom <- Racmacs:::mapPlotLims(map, sera = TRUE)
xlim_no_zoom <- round(lims_no_zoom$xlim)
ylim_no_zoom <- round(lims_no_zoom$ylim)
lims_zoom <- Racmacs:::mapPlotLims(map, sera = FALSE)
xlim_zoom <- round(lims_zoom$xlim)
ylim_zoom <- round(lims_zoom$ylim)
# plot map
layout(matrix(c(1:2), ncol = 2, byrow = FALSE))
par(oma=c(0, 0, 0, 0), mar=c(0, 0.5, 0, 0))
plot(map, plot_labels = "antigens", label.offset = 1.5, xlim = xlim_zoom, ylim = ylim_zoom, plot_stress = TRUE)
plot(map, plot_labels = "antigens", label.offset = 1.5, xlim = xlim_no_zoom, ylim = ylim_no_zoom, plot_stress = TRUE)
Above we see the same antigenic map, once zoomed in (left) and once zoomed out (right). Variants are shown as coloured circles, sera are shown as open squares in the colour of their root variant (or as triangles pointing toward their position if they outside of the shown map area). Characteristics of a single exposure map are that sera are located close to their infecting/homologous variant as they have the highest titer against this variant. Hence, the antigenic distance to this variant is 0. All other variants are located at the Euclidean distance that best matches the log2(fold change) of titers from the homologous variant to this variant. However, you will see that sera are not at exactly the same position as their homologous variant. That is because the optimization algorithm finds the optimum for all serum and variant coordinates to best match the titer fold changes. As there is titer variation between individual samples, there will be variation in their positions. The mismatch between target titer distance (titer fold changes) and euclidean map distance (antigenic distance) is given by the map stress in the lower left corner. The lower the stress, the better the match between target and map distance.
Racmacs provides an interactive viewer with featuers such as map randomization, reoptimization, movement of individual objects, map stress investigation. To properly use these featuers, non-positioned sera and variants (due to too many not measured/undetectable) titers need to be removed.
If the viewer does not properly render here, you can open it in any R session.
remove_na_sera <- function(map){
na_sera <- srCoords(map)
na_sera <- rownames(na_sera)[is.na(na_sera)[,1]]
map <- removeSera(map, na_sera)
return(map)
}
remove_na_antigens <- function(map){
na_sera <- agCoords(map)
na_sera <- rownames(na_sera)[is.na(na_sera)[,1]]
map <- removeAntigens(map, na_sera)
return(map)
}
remove_na_coords <- function(map){
map <- remove_na_sera(map)
map <- remove_na_antigens(map)
return(map)
}
Racmacs::view(remove_na_coords(map))
The plot below shows how much each variant and serum is allowed to move without increasing the map stress by more than one unit. It illustrates that more >LOD titers in different types of sera increase the positional resolution (e.g. the blue D614G has >LOD in almost all sera, whereas BA.1 has <LOD in almost all sera). For this diagnostic step, again sera that could not be positioned need to be removed from the map.
# triangulation does not work with not-positioned sera - sera cannot be position if too many titers are <LOD
plot(triangulationBlobs(relaxMap(remove_na_coords(map)), stress_lim = 1, grid_spacing = 0.05), fill.alpha = 0.9, plot_labels = "antigens", outline.alpha = 0.9, xlim = xlim_zoom, ylim = ylim_zoom,
grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
To demonstrate the importance of different types of sera on map conformation and positional resolution, I will now create maps subset to
only 2x Pfizer vaccine sera,
previous + alpha sera (different homologous variant but similar antigenic profile)
previous + beta sera (sera with more distinct antigenic profile)
previosu + BA.1 sera (sera with very distinct antigenic profile)
subset_and_optimize_map <- function(full_map, kept_sr_groups){
map <- removeSera(full_map, srNames(full_map)[!as.character(srGroups(full_map)) %in% kept_sr_groups])
map <- optimizeMap(
map,
number_of_dimensions = 2,
number_of_optimizations = 1000,
minimum_column_basis = "none",
options = list(ignore_disconnected = TRUE)
)
map <- realignMap(relaxMap(remove_na_coords(map)), full_map)
srOutlineWidth(map) <- 2
srSize(map) <- 9*4
agSize(map) <- 18*4
srOutline(map) <- adjustcolor(srOutline(map), alpha.f = 0.7)
return(map)
}
vax_map <- subset_and_optimize_map(map, c("BNT/BNT"))
alpha_map <- subset_and_optimize_map(map, c("BNT/BNT", "alpha conv."))
beta_map <- subset_and_optimize_map(map, c("BNT/BNT", "alpha conv.", "beta conv."))
ba1_map <- subset_and_optimize_map(map, c("BNT/BNT", "alpha conv.", "BA.1 conv.", "beta conv."))
ylim_no_zoom[2] <- ylim_no_zoom[2] + 1
layout(matrix(c(1:8), ncol = 4, byrow = FALSE))
par(oma=c(0, 0, 0, 0), mar=c(0.5, 0, 0.1, 0))
plot(procrustesMap(vax_map, map, sera = FALSE), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_no_zoom, ylim = ylim_no_zoom, grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(triangulationBlobs(vax_map, stress_lim = 1, grid_spacing = 0.05), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_no_zoom, ylim = ylim_no_zoom,
grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(procrustesMap(alpha_map, map, sera = FALSE), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_no_zoom, ylim = ylim_no_zoom, grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(triangulationBlobs(alpha_map, stress_lim = 1, grid_spacing = 0.05), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_no_zoom, ylim = ylim_no_zoom,
grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(procrustesMap(beta_map, map, sera = FALSE), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_no_zoom, ylim = ylim_no_zoom, grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(triangulationBlobs(beta_map, stress_lim = 1, grid_spacing = 0.05), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_no_zoom, ylim = ylim_no_zoom,
grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(procrustesMap(ba1_map, map, sera = FALSE), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_no_zoom, ylim = ylim_no_zoom, grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(triangulationBlobs(ba1_map, stress_lim = 1, grid_spacing = 0.05), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_no_zoom, ylim = ylim_no_zoom,
grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
From left to right:
only 2x Pfizer vaccine sera,
previous + alpha sera (different homologous variant but similar antigenic profile)
previous + beta sera (sera with more distinct antigenic profile)
previosu + BA.1 sera (sera with very distinct antigenic profile)
The top row shows the antigenic maps with arrows pointing to the variants’ positions in the full map, the bottom row shows the constant force loci (triangulation blobs).
The third map gives a very similar confirmation as the “full map” but the positional resolution of the Omicron variants is very low.
Multiple exposures increase cross-neutralization and hence fold changes of titers between variants are less pronounced. Creating a map from such sera results in a map where variants are closer together than in single exposure sera maps and an underestimation of antigenic relationships. Multiple exposures do not change the antigenic relationships of virus strains, but change how cross-reactive sera are. We can see that by comparing fold changes from the homologous titer variant after 2 and 3 Pfizer vaccines:
# get GMT fold change from D614G
gmt_data %>%
filter(sr_group %in% c("BNT/BNT", "BNT/BNT/BNT")) %>%
group_by(sr_group) %>%
mutate(log_fc = logtiter[ag_name == "D614G"] - logtiter,
fc = ifelse(log_fc > 0, paste0("-", round(2^abs(log_fc), 1)), round(2^abs(log_fc), 1))) -> fc_data
do_titerplot(data_long,
target_sr_groups = c("BNT/BNT", "BNT/BNT/BNT"),
gmt_data,
sr_group_colors,
ag_order = rownames(titer_data)) +
geom_text(data = fc_data, aes(label = fc), y = 10.5) +
geom_text(data = fc_data, aes(label = abs(round(log_fc, 1))), y = 10)
The top row of numbers lists the GMT fold change from D614G, the row below the log2(fold change) and thereby the antigenic distance. The antigenic distance difference is the most pronounced for the most distant Omicron variants.
To illustrate this effect on an antigenic map, we will create a map from 3x Pfizer vaccine sera and compare it to the 2x Pfizer vaccine sera:
# make table for 3x Pfizer sera
multi_exposure_data %>%
filter(sr_group == "BNT/BNT/BNT") %>%
select(!sr_group) %>%
pivot_wider(names_from = "sr_name", values_from = "titer") %>%
column_to_rownames("ag_name") -> table_3x
# make 3x map
# create the acmap object
map_3x <- acmap(
ag_names = rownames(table_3x),
sr_names = colnames(table_3x),
titer_table = table_3x)
# The dilution stepsize gives the dilution steps of the neutralization assay on the log2 scale, e.g.: 0 = continuous read out, 1 = two-fold dilutions
dilutionStepsize(map_3x) <- 0
# optimize the map in 2 dimensions
map_3x <- optimizeMap(
map_3x,
number_of_dimensions = 2,
number_of_optimizations = 1000,
minimum_column_basis = "none",
options = list(ignore_disconnected = TRUE) # Map contains disconnected points (points that are not connected through any path of detectable titers so cannot be coordinated relative to each other)
)
map_3x <- realignMap(map_3x, map)
srOutline(map_3x) <- "grey20"
agFill(map_3x) <- mapColors[agNames(map_3x), "Color"]
srOutlineWidth(map_3x) <- 2
srOutline(map_3x) <- adjustcolor(srOutline(map_3x), alpha.f = 0.7)
agSize(map_3x) <- 18*4
srSize(map_3x) <- 9*4
# Do plotting
layout(matrix(c(1:3), ncol = 3, byrow = FALSE))
par(oma=c(0, 0, 0, 0), mar=c(0.1, 0.5,0, 0))
plot(map_3x, fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_zoom, ylim = ylim_zoom + 0.5, grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(realignMap(vax_map, map_3x), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_zoom, ylim = ylim_zoom+0.5, grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
plot(procrustesMap(map_3x, vax_map, sera = FALSE), fill.alpha = 0.9, plot_labels = FALSE, outline.alpha = 0.9, xlim = xlim_zoom, ylim = ylim_zoom + 0.5, grid.col = "#cfcfcf", grid.margin.col="#7d7d7d", cex=0.3, plot_stress = TRUE)
On the left, we see the map from 3x Pfizer vaccine sera, the middle map is the 2x Pfizer vaccine map and on the right, the 3x Pfizer map is shown with arrows pointing to the variants’ positions in the 2x Pfizer map. The arrows point outside, indicating that the variants are closer to each other in the 3x Pfizer map, and therefore look more antigenically similar than in the 2x Pfizer map. However, the virus isolates used in the assay are the same. Again, it is not the antigenic relationships and properties of the viruses that change with more exposures, but the antibody composition in the multi-exposure sera.
Antibody landscapes better represent multiple exposure sera. Neutralization titers are mapped as a continuous surface above an antigenic map: Titers against the respective variant are mapped above the variant in a third dimension, the height representing the titer magnitude, and connected through a surface, the slope of which indicates the cross-reactivity. The smaller the slope, the more cross-reactive the sample.
In the single cone landscape approached used here, cone coordinates (x, y) and cone heights (z) are fitted to titers per subject. The slope of this cone is optimised over all subjects per serum group, assuming the same decrease of neutralization across antigenic space. The cone slopes indicate the breadth of neutralization titers. As immune history becomes more complex, for example for influenza, a single cone approach does not provide a good approach and a loess fit should be used.
The description above focusses on antibody neutralization, not binding. The landscape fitting approach remains the same, however.
# orientation for landscapes
angle <- list(
rotation = c(-1.4681,0.004,-0.0162),
translation = c(0, 0.05,0.1),
zoom = 1.45
)
# a function to get titer tables from long format to matrix form as input for the antibody landscape fit
get_titertable <- function(data, group) {
data %>%
select(
ag_name,
sr_name,
titer
) %>%
mutate(
titer = replace(titer, is.na(titer), "*")
) %>%
pivot_wider(
id_cols = sr_name,
names_from = ag_name,
values_from = titer
) %>%
as.matrix() -> titermatrix
attr(titermatrix, "sr_group") <- group$sr_group
rownames(titermatrix) <- titermatrix[,"sr_name"]
titermatrix <- titermatrix[,-1]
return(titermatrix)
}
# plot the base map as r3js object
base_plot_data3js <- function(map, lndscp_fits, highlighted_ags, lims, ag_plot_names,
add_border = TRUE, add_axis = TRUE, add_ag_labels = TRUE){
# get coordinates for variants that should be plotted (highlighted ags)
x_coords <- c(agCoords(map)[agNames(map) %in% highlighted_ags, 1])
y_coords <- c(agCoords(map)[agNames(map) %in% highlighted_ags, 2])
z_coords <- rep(0.02, length(highlighted_ags))
ag_point_size <- c(rep(14, length(highlighted_ags))) / 5
ag_col <- c(agOutline(map)[agNames(map) %in% highlighted_ags])
ag_fill <- c(agFill(map)[agNames(map) %in% highlighted_ags])
labels <- c(ag_plot_names[agNames(map) %in% highlighted_ags])
border_col <- "grey50"
z_lims <- c(0,10)
axis_at <- seq(z_lims[1], z_lims[2],2)
# Setup plot
data3js <- ablandscapes:::lndscp3d_setup(
xlim = lims$xlim,
ylim = lims$ylim,
zlim = z_lims,
aspect.z = 0.5,
options = list(
lwd.grid = 0.05,
sidegrid.lwd = 1,
sidegrid.col = border_col,
sidegrid.at = list("z" = axis_at),
zaxt = "log"
),
show.axis = FALSE
)
# add z axis
if(add_axis){
axis_labels <- 2^axis_at*10
data3js <- r3js::axis3js(
data3js,
side = "z",
at = axis_at,
labels = axis_labels,
cornerside = "f",
size = 20,
alignment = "right"
)
}
# Add basemap
data3js <- lndscp3d_map(
data3js = data3js,
fit = lndscp_fits[[1]],
xlim = lims$xlim,
ylim = lims$ylim,
zlim = c(0, 10),
show.map.sera = FALSE,
options = list(
opacity.basemap = 0.3
)
)
# add variants
data3js <- r3js::points3js(
data3js,
x = x_coords,
y = y_coords,
z = z_coords,
size = ag_point_size,
col = ag_col,
fill = ag_fill,
lwd = 0.5,
opacity = 1,
highlight = list(col = "red"),
label = labels,
toggle = "Basepoints",
depthWrite = FALSE,
shape = "circle filled"
)
if(add_ag_labels){
text_x <- x_coords
text_y <- c(y_coords - ag_point_size*0.2)
data3js <- r3js::text3js(
data3js,
x = text_x,
y = text_y,
z = z_coords,
text = labels,
size = c(rep(10*0.02, length(text_x))),
alignment = "center"
)
}
# thicker border of the base map
if(add_border){
data3js <- lines3js(data3js, x = c(lims$xlim[1],lims$xlim[1]), y = c(lims$ylim[1], lims$ylim[2]), z = c(0, 0),
lwd = 1.2, col = border_col)
data3js <- lines3js(data3js, x = c(lims$xlim[2],lims$xlim[2]), y = c(lims$ylim[1], lims$ylim[2]), z = c(0, 0),
lwd = 1.2, col = border_col)
# y border
data3js <- lines3js(data3js, x = c(lims$xlim[1],lims$xlim[2]), y = c(lims$ylim[1], lims$ylim[1]), z = c(0, 0),
lwd = 1.2, col = border_col)
data3js <- lines3js(data3js, x = c(lims$xlim[1],lims$xlim[2]), y = c(lims$ylim[2], lims$ylim[2]), z = c(0, 0),
lwd = 1.2, col = border_col)
data3js <- r3js::box3js(
data3js,
col = border_col
)
}
return(data3js)
}
# plot landscapes from a list that contains the fits
plot_landscapes_from_list <- function(data3js, # base map as r3js object
titertables_groups, # grouped titer table object for serum groups
lndscp_fits, # landscape fits as list
map, # map for x,y coordinates of GMT points
highlighted_ags, # which map antigens should be highlighted (GMT coordinates)
lndscp_colors, # colors for the landscapes
gmt_data = NULL,
show_gmts = TRUE,
show.individual.surfaces = FALSE,
options.individual.surfaces = list(opacity.surface.grid = 0.4,
opacity.surface = 0.2,
col.surface = "grey70",
col.surface.grid = "grey70"),
show_landscapes = TRUE){
# get x and y coords for gmt points
x_coords <- c(agCoords(map)[agNames(map) %in% highlighted_ags, 1])
y_coords <- c(agCoords(map)[agNames(map) %in% highlighted_ags, 2])
coords <- cbind(x_coords, y_coords)
coords <- coords[!is.na(x_coords),]
# for each landscape fit object (serum group), do the plotting
for (i in seq_along(lndscp_fits)) {
# select sr group
srg <- titertables_groups$sr_group[i]
lndscp_fit <- lndscp_fits[[i]]
for (j in seq_len(nrow(coords))) {
if(show_gmts){
if(is.null(gmt_data)){
warning("Error: You want to plot GMT but did not provide GMT data")
return()
}
# select the gmt data for the serum group
gmts <- gmt_data %>%
filter(sr_group == srg)
gmts <- gmts[match(rownames(coords), gmts$ag_name),]
# plot GMTs
data3js <- r3js::lines3js(
data3js,
x = rep(coords[j, 1], 2),
y = rep(coords[j, 2], 2),
z = c(0, gmts$gmt[j]),
col = "grey50",
toggle = sprintf("GMT, %s", srg),
geometry = TRUE,
opacity = 0.7,
lwd = 0.2
)
data3js <- r3js::points3js(
data3js,
x = coords[j, 1],
y = coords[j, 2],
z = gmts$gmt[j],
size = 0.9,
col = lndscp_colors[srg, "Color"],
toggle = sprintf("GMT, %s", srg),
opacity = 1
)
}
}
if(show_landscapes){
# Add GMT landscapes
data3js <- lndscp3d_surface(
data3js = data3js,
object = lndscp_fit,
crop2chull = FALSE,
toggle = sprintf("GMT landscape, %s", srg),
grid_spacing = 0.5,
padding = 0.2,
options = list(
col.surface = lndscp_colors[srg, "Color"],
opacity.surface = 1
)
)
fit <- lndscp_fit
# show individual landscapes
if (show.individual.surfaces) {
for (i in seq_len(nrow(fit$titers))) {
individual_fit <- fit
individual_fit$titers <- fit$titers[i, ]
individual_fit$logtiters <- fit$logtiters[i, ]
individual_fit$logtiters.upper <- individual_fit$logtiters.upper[i,
]
individual_fit$logtiters.lower <- individual_fit$logtiters.lower[i,
]
individual_fit$lessthans <- individual_fit$lessthans[i,
]
individual_fit$morethans <- individual_fit$morethans[i,
]
individual_fit$fitted.values <- NULL
individual_fit$residuals <- NULL
individual_fit$residuals.lessthan <- NULL
individual_fit$residuals.morethan <- NULL
if (!is.null(individual_fit$cone)) {
individual_fit$cone$cone_coords <- individual_fit$cone$cone_coords[i,
, drop = F]
individual_fit$cone$cone_heights <- individual_fit$cone$cone_heights[i]
}
data3js <- lndscp3d_surface(data3js = data3js, object = individual_fit,
crop2chull = FALSE, grid_spacing = 0.5,
padding = 0.1,
options = options.individual.surfaces,
toggle = "Individual landscapes")
}
}
}
}
return(data3js)
}
# set orientation
set_r3js_orentation <- function(data3js, angle){
r3js(
data3js,
rotation = angle$rotation,
zoom = angle$zoom
)
}
# group multi exposure data by serum group
multi_exposure_data %>%
group_by(
sr_group
) -> titerdata
titerdata %>%
group_map(
get_titertable
) -> titertables
# do the fit per serum group on base map
lndscp_fits <- lapply(titertables, function(titer_table){
ablandscape.fit(
titers = titer_table, #errors can occur if the titertable for one serum group contains only one sample (then this sample is passed as a vector rather than a matrix) and if all titers in a sample are <LOD/NA. The titers are fitted, not the logtiters. <LOD values are estimated, but if all values are <LOD there is too little data to estimate from
bandwidth = 1,
degree = 1,
method = "cone",
error.sd = 1,
acmap = map,
control = list(
optimise.cone.slope = TRUE
)
)
})
# calculate GMTs per serum group
titertables_groups <- group_data(titerdata)
# plot the base map
data3js <- base_plot_data3js(map, lndscp_fits, highlighted_ags = agNames(map), ag_plot_names = agNames(map), lims = lims_zoom)
set_r3js_orentation(data3js, angle)
# plot 2x vax gmts
pfizer2_ind <- match("BNT/BNT", titertables_groups$sr_group)
# do lndscp colors in correct format
lndscp_colors <- sr_group_colors %>%
column_to_rownames("SerumGroup")
gmt_2x <- plot_landscapes_from_list(data3js = data3js, titertables_groups = titertables_groups[pfizer2_ind,], lndscp_fits = lndscp_fits[pfizer2_ind], map = map, gmt_data = gmt_data, highlighted_ags = agNames(map), lndscp_colors = lndscp_colors, show_gmts = TRUE, show.individual.surfaces = FALSE, show_landscapes = FALSE)
set_r3js_orentation(gmt_2x, angle)
gmt_scape_2x <- plot_landscapes_from_list(data3js = data3js, titertables_groups = titertables_groups[pfizer2_ind,], lndscp_fits = lndscp_fits[pfizer2_ind], map = map, gmt_data = gmt_data, highlighted_ags = agNames(map), lndscp_colors = lndscp_colors, show_gmts = TRUE, show.individual.surfaces = FALSE, show_landscapes = TRUE)
set_r3js_orentation(gmt_scape_2x, angle)
gmt_scape_2x <- plot_landscapes_from_list(data3js = data3js, titertables_groups = titertables_groups[pfizer2_ind,], lndscp_fits = lndscp_fits[pfizer2_ind], map = map, gmt_data = gmt_data, highlighted_ags = agNames(map), lndscp_colors = lndscp_colors, show_gmts = TRUE, show.individual.surfaces = TRUE, show_landscapes = TRUE)
set_r3js_orentation(gmt_scape_2x, angle)
target_groups <- grep("BNT/BNT|BNT/BNT/BNT", titertables_groups$sr_group)[c(1,3)]
lndscps_plot <- plot_landscapes_from_list(data3js = data3js, titertables_groups = titertables_groups[target_groups,], lndscp_fits = lndscp_fits[target_groups], map = map, gmt_data = gmt_data, highlighted_ags = agNames(map), lndscp_colors = lndscp_colors, show_gmts = TRUE, show_landscapes = TRUE)
set_r3js_orentation(lndscps_plot, angle)
The antibody landscapes after different doses illustrate the
increasing breadth after repeated exposure. The cone slope of each
landscape is a measure of its breadth and can be accessed via the
cone field of each antibody landscape fit, which also
contains the fitted cone coordinates and cone heights per subject. The
underlying geometry of the antigenic map - one antigenic distance unit
corresponds to one two-fold change of neutralization titers - means that
single exposure sera, from which the map was built, have a landscape
slope around 1 on the log2 scale: With each antigenic distance unit away
from the cone apex, the titers decrease by one two-fold. The slopes of
multi-exposure sera are much lower than the slopes of single exposure
sera:
slopes <- lapply(lndscp_fits, function(fit){
round(fit$cone$cone_slope,1)
})
slopes_df <- data.frame("Serum group" = factor(titertables_groups$sr_group, levels = levels(srGroups(map))),
"Landscape slope" = unlist(slopes)) %>%
arrange(Serum.group)
kable(slopes_df, col.names = c("Serum group", "Landscape slope")) %>%
kable_styling(full_width = F)
| Serum group | Landscape slope |
|---|---|
| mRNA1273/mRNA1273 | 1.2 |
| BNT/BNT | 1.2 |
| AZ/BNT | 1.1 |
| AZ/AZ | 1.1 |
| WT conv. | 1.0 |
| alpha conv. | 1.3 |
| beta conv. | 1.5 |
| delta conv. | 1.2 |
| BA.1 conv. | 0.8 |
| BA.2 conv. | 1.3 |
| Vacc+BA.1 | 0.7 |
| Vacc+BA.2 | 0.5 |
| BA.1 reinf. | 0.6 |
| BA.2 reinf. | 0.4 |
| BNT/BNT/BNT | 0.7 |
| Vacc+BA.1 reinf. | 0.7 |
| AZ/AZ+delta | 0.8 |
| BNT/BNT+delta | 0.8 |
lndscps_plot <- plot_landscapes_from_list(data3js = data3js, titertables_groups = titertables_groups, lndscp_fits = lndscp_fits, map = map, gmt_data = gmt_data, highlighted_ags = agNames(map), lndscp_colors = lndscp_colors, show_gmts = TRUE, show_landscapes = TRUE)
set_r3js_orentation(lndscps_plot, angle)
An easy first visual inspection is also the direct comparison between landscape fitted GMT and measured GMT:
# bind landscape fit and calculated gmt
combine_landscape_and_calculated_gmt <- function(lndscp_fits, gmt_data, sr_group_fields = 2){
lndscp_gmts <- lapply(lndscp_fits, function(x){
gmts <- data.frame(logtiter = x$fitted.value,
ag_name = names(x$fitted.value),
sr_group = str_split(rownames(x$logtiters)[1], "_")[[1]][sr_group_fields])
return(gmts)
})
lndscp_gmts <- do.call(rbind, lndscp_gmts)
## compare lndscp gmts and claculated gmts
comb_gmt <- rbind(lndscp_gmts %>%
mutate(Data = "Fitted Landscape GMT"),
gmt_data %>%
select(sr_group, ag_name, logtiter) %>%
unique() %>%
mutate(Data = "Calculated GMT"))
return(comb_gmt)
}
comb_data <- combine_landscape_and_calculated_gmt(lndscp_fits, gmt_data %>%
mutate(logtiter = gmt) %>%
select(!gmt),
sr_group_fields = 2)
# plot comparison
comb_data %>%
mutate(sr_group = factor(sr_group, levels = rownames(lndscp_colors))) %>%
ggplot(aes(x = ag_name, y = logtiter, color = Data, fill = Data)) +
geom_line(aes(group = Data), position = position_dodge(width = 0.05)) +
geom_point(shape = 21, color = "grey20", position = position_dodge(width = 0.05)) +
scale_x_discrete(name = "Variant",
limits = agNames(map)) +
scale_y_continuous(breaks = seq(log2(16/10), 12, by = 2),
labels = function(x) ifelse(x == log2(16/10), paste0("<", 16), round(2^x*10)),
limits = c(-3, 12),
name = "GMT") +
annotate(
"rect",
xmin= -Inf,
xmax = Inf,
ymin = -Inf,
ymax = log2(16/10),
fill = "grey50",
color = NA,
alpha = 0.2
) +
facet_wrap(~sr_group,
labeller = label_wrap_gen(multi_line = TRUE)) +
theme_bw() +
theme(strip.background = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
axis.text.y = element_text(size = 8)) -> p
p
Each ablandscape fit object contains a residuals field
to check for error between measured logtiter and fitted landscape GMT
for individual samples. The residuals are given on the log2 scale, so a
residual of 1 corresponds to 1 2-fold change difference between measured
titer and fitted landscape GMT.
# get residuals into long format
residuals_to_long <- function(residuals, values_name = "residuals"){
ag_names <- colnames(residuals)
as.data.frame(residuals) %>%
rownames_to_column(var = "sr_name") %>%
pivot_longer(cols = all_of(ag_names), names_to = "ag_name", values_to = values_name) -> residuals_long
return(residuals_long)
}
# extract residuals from landscape fit function
combine_residuals <- function(fit, sr_group_fields = 2){
residuals <- fit$residuals
less_thans <- fit$residuals.lessthan
more_thans <- fit$residuals.morethan
long_res <- residuals_to_long(residuals, "residuals")
long_less <- residuals_to_long(less_thans, "less_than")
long_more <- residuals_to_long(more_thans, "more_than")
serum_group <- paste0(str_split(long_res$sr_name[1], "_")[[1]][sr_group_fields], collapse = "_")
comb <- long_res %>%
left_join(long_less, by = c("ag_name", "sr_name")) %>%
left_join(long_more, by = c("ag_name", "sr_name")) %>%
mutate(measured = ifelse(is.na(residuals), ifelse(is.na(less_than), "more_than", "less_than"), "detectable"),
residuals = ifelse(is.na(residuals), ifelse(is.na(less_than), more_than, less_than), residuals)) %>%
select(!less_than:more_than) %>%
mutate(sr_group = serum_group)
return(comb)
}
# get root mean square error per variant
rmse_per_variant <- function(lndscp_fits, sr_group_fields = 2){
all_residuals <- lapply(lndscp_fits, function(x) combine_residuals(x, sr_group_fields))
all_residuals <- do.call(rbind, all_residuals)
all_residuals %>%
filter(!is.na(residuals)) %>%
group_by(sr_group) %>%
summarize(ag_name = "Total",
rmse = sqrt(sum(residuals^2, na.rm = TRUE)/(length(residuals))), # Residual standard error
residual_type = "All variants") -> total_ag
all_residuals %>%
filter(!is.na(residuals)) %>%
group_by(ag_name, sr_group) %>%
mutate(rmse = sqrt(sum(residuals^2, na.rm = TRUE)/(length(residuals))),
residual_type = "By variant") %>%
plyr::rbind.fill(., total_ag) -> ssr
return(ssr)
}
residuals <- rmse_per_variant(lndscp_fits, sr_group_fields = 2) %>%
mutate(sr_group = factor(sr_group, levels = levels(srGroups(map))))
# plot residuals
residuals %>%
ggplot(aes(x = ag_name, y = rmse, fill = ag_name)) +
geom_point(color = "grey20", shape = 21) +
scale_x_discrete(limits = agNames(map),
name = "Variant") +
ylab("RMSE (log2(FC))") +
ylim(c(0,max(ceiling(residuals$rmse)))) +
facet_wrap(~sr_group) +
scale_fill_manual(values = setNames(mapColors$Color, rownames(mapColors)),
name = "Variant") +
theme_bw() +
theme(strip.background = element_blank(),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)) -> rmse_plot_variant
# plot residuals
residuals %>%
filter(ag_name == "Total") %>%
ggplot(aes(x = sr_group, y = rmse, fill = sr_group)) +
geom_point(color = "grey20", shape = 21) +
ylab("RMSE (log2(FC))") +
ylim(c(0,max(ceiling(residuals$rmse)))) +
scale_fill_manual(values = setNames(mapColors$Color, rownames(mapColors)),
name = "Serum group") +
theme_bw() +
xlab("Serum group") +
theme(strip.background = element_blank(),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)) -> rmse_sr_groups
rmse_sr_groups + rmse_plot_variant + plot_layout(widths = c(1, 1.75))
Sometimes virus variants can exhibit high reactivity in an assay, meaning that measured titers against this variant are high due to assay, but not antigenic effects. It can also be the other way around, that a variant is low reactivity and titers against this variant are low in all samples. Judging whether a variant exhibits a reactivity bias is not always straight forward and requires external information, such as sequence substitutions.
The data set above contains the Gamma (P.1.1) variant, which was excluded in the demonstration above, because it exhibits a reactivity bias. Two factors in the current data set point towards this: Firstly, Gamma titers are higher than homologous titers in many sera, also in sera for which the homologous variant has distinct amino acids at antigenically important sites. Secondly, titers against the Beta (B.1.351) variant, which differs in spike from Gamma at position 417, are much lower in the same sera:
# read in data for multiple exposure landscapes
multi_exposure_data <- multi_exposure_data_b
multi_exposure_data %>%
# log2 transform the data and set <LOD values to LOD/2
mutate(logtiter = ifelse(titer == "<16", log2(16/20), log2(as.numeric(titer)/10)),
sr_group = factor(sr_group, levels = sr_group_colors$SerumGroup)) -> data_long
# calculate
data_long %>%
group_by(
sr_group,
ag_name
) %>%
# if you want to use quick LOD/2 method for GMT calculation
summarize(logtiter = mean(logtiter, na.rm = TRUE)) %>%
mutate(logtiter = ifelse(logtiter < log2(0.8), log2(0.8), logtiter),
sr_name = "GMT",
gmt = logtiter)-> gmt_data
# plot the data
do_titerplot(data_long, target_sr_groups = unique(data_long$sr_group),
gmt_data = gmt_data,
sr_group_colors = sr_group_colors,
ag_order = agNames(og_map)) +
annotate("rect",
xmin = "P.1.1",
xmax = "B.1.351",
ymin = -Inf,
ymax = log2(0.8),
fill = "red")
The Racmacs package has a function to deal with antigen reactivity biases. Reactivity adjustment are passed as a named vector on the log2 scale. So a reactivity adjustment of -1 means that all titers for this variant are multiplied with a factor of 0.5 (2^-1).
map <- og_map
# reduce P.1.1 reactivity by 1 dilution step. Alternatively, you can optimize an antigens reactivity using the function below. 0 = no adjustment, NA = optimization, value for numeric adjustment on the log2 scale
ag_reactivities <- rep(0, length(agNames(map)))
ag_reactivities[agNames(map) == "P.1.1"] <- -1
map <- optimizeAgReactivity(map, fixed_ag_reactivities = ag_reactivities, reoptimize = TRUE)
The adjusted titer data can be accessed using Racmacs’
adjustedTiterTable() and
adjustedLogTiterTable() functions.